%title:             insurance model with misperceived health type correlation
%                   with actual cost = 0
%author:            neale mahoney
%lastest revision:  6/29/14


%% Setup
clear all; clc;

% cd(<folder directory>);

version = 'v2';

%% Primatives

fprintf('Displaying primatives \n\n');

n_i = 16e3;             %number of individuals
n_c = 1e3;              %number of cost shocks per individual
smooth_type = 'lowess'; %moving, lowess, rlowess
bw = .1;                %bandsmooth for moving average smoother

%% Joint distribution of risk aversion and expected medical costs

%distribution of absoluate risk aversion from Table 3 of Handel, Hendel and Whinston (2013)
m = 4.39e-4;
v = (6.63e-5)^2;
mu_alpha = log((m^2)/sqrt(v+m^2));
sigma_alpha = sqrt(log(v/(m^2)+1));

%gamble interpretation
f = zeros(100,1);
for i = 1:100,
    f(i) = cara_function(0,m,[-100 i],[0 0]);
end

[value index] = min(abs(f));
fprintf('Mean risk adversion makes individual different between 50-50 gamble for (100,%i) and 0 \n', index)
    
% distribution of health type
m = .9794951;
v = 1.377593^2;
mu_l = log((m^2)/sqrt(v+m^2));
sigma_l = sqrt(log(v/(m^2)+1));

%joint normal of risk preference, actual health type, and perceived health type
mu = [mu_alpha; mu_l; mu_l];

%correlation between actual and perceived health type
rho = 0;
cov = rho*sigma_l*sigma_l;

%assume risk aversion and actual/perceived health type are uncorrelated in the population
sigma = [sigma_alpha^2 0 0; 0 sigma_l^2 cov; 0 cov sigma_l^2];
X = mvnrnd(mu,sigma,n_i);
 
%distribution of risk preferences and health type
alpha = exp(X(:,1));
lambda = exp(X(:,2));

%distribution of realized costs
m = 3138.765; 
v = 10125.74^2;
mu_c = log((m^2)/sqrt(v+m^2));
sigma_c = sqrt(log(v/(m^2)+1));
rho = .498481; 

%distribution of realized cost conditional on health type
c = zeros(n_i,n_c);
part_1 = mu_c;
part_2 = rho*sigma_c*(X(:,2)-mu_l)/sigma_l;

for i = 1:n_i   
    part_3 = sqrt((1-rho^2))*sigma_c*randn(1,n_c);
    c(i,:) = exp(part_1 + part_2(i) + part_3);
end

%distribution of perceived costs
c_per = zeros(n_i,n_c);
part_1 = mu_c;
part_2 = rho*sigma_c*(X(:,3)-mu_l)/sigma_l;

for i = 1:n_i   
    part_3 = sqrt((1-rho^2))*sigma_c*randn(1,n_c);
    c_per(i,:) = exp(part_1 + part_2(i) + part_3);
end

%display data
fprintf('Mean absolute risk aversion parameter is %f \n',mean(alpha));
fprintf('Std dev of absolute risk aversion parameter is %f \n',std(alpha));

fprintf('Mean population cost is %i \n',round(mean(reshape(c,n_i*n_c,1))));
fprintf('Std dev of cost is %i \n',round(std(reshape(c,n_i*n_c,1))));

fprintf('Mean perceived cost is %i \n',round(mean(reshape(c_per,n_i*n_c,1))));
fprintf('Std dev of cost is %i \n',round(std(reshape(c_per,n_i*n_c,1))));

c_mean = mean(c,2);
c_per_mean = mean(c_per,2);
fprintf('\nCorrelation between acutal and perceived mean costs is %i \n', corr(c_mean, c_per_mean));
%scatter(c_mean, c_per_mean);

%av 60 plan: low quality plan that covers 60% of medical cost on average
c_oop_60 = min(min(c,500+.4*(c-500)),8e3);  %500 deductible, 8000 out of pocket maximum
c_oop_60_per = min(min(c_per,500+.4*(c_per-500)),8e3);

c_ins_60 = c - c_oop_60;
c_ins_mean_60 = mean(c_ins_60,2);

fprintf('AV 60 mean plan cost is %i \n',round(mean(c_ins_mean_60)))
fprintf('AV 60 mean oop cost is %i \n',round(mean(reshape(c_oop_60,n_i*n_c,1))))
fprintf('AV 60 AV is %f \n',mean(c_ins_mean_60)/mean(reshape(c,n_i*n_c,1)))

% av 90 plan: high quality plan that covers 90% of medical cost on average
c_oop_90 = min(.1*c,8e3);   %no deductible, 8000 out of pocket maximum
c_oop_90_per = min(.1*c_per,8e3);

c_ins_90 = c - c_oop_90;
c_ins_mean_90 = mean(c_ins_90,2);

fprintf('AV 90 mean plan cost is %i \n',round(mean(c_ins_mean_90)))
fprintf('AV 90 mean oop cost is %i \n',round(mean(reshape(c_oop_90,n_i*n_c,1))))
fprintf('AV 90 AV is %f \n',mean(c_ins_mean_90)/mean(reshape(c,n_i*n_c,1)))

%allow bankruptcy to cap out of pocket risk
c = min(c,20e3);
c_per = min(c_per,20e3);

%% wtp
fprintf('\n\nEstimating WTP \n\n');

v_grid = 0:10:5e4;
v = zeros(n_i,1);
for i = 1:n_i,
    %v: willingness to pay additional premium such that 
    %the customer is indifferent between actual costs of high and low quality plan
    if mod(i,25) == 0, fprintf('Estimating WTP for observation %i \n',i); end
    v(i) = wtp_function(alpha(i),v_grid,c_oop_90(i,:),c_oop_60(i,:));

end

%% wtp with noise (Handel, Kolstad, and Spinjiwin)
fprintf('\n\nEstimating WTP \n\n');

v_grid = 0:10:5e4;
v_per = zeros(n_i,1);
for i = 1:n_i,
    %v: willingness to pay additional premium such that 
    %the customer is indifferent between perceived costs of high and low quality plan
    if mod(i,25) == 0, fprintf('Estimating perceived WTP for observation %i \n',i); end
    v_per(i) = wtp_function(alpha(i),v_grid,c_oop_90_per(i,:),c_oop_60_per(i,:));

end


%% Equilibrium - Baseline 

%sort matrix X in descending order by wtp and smooth
wtp_shift = 900;
X = [v + wtp_shift c_ins_mean_90-c_ins_mean_60];
X = sortrows(X,-1);

v_sort = smooth(X(:,1),bw,smooth_type);
mc = smooth(X(:,2),bw,smooth_type);
q = transpose((1:n_i)./n_i);

plot_switch = 1;
plot_number = 2;
plot_dir = '';
plot_name = 'all_rho_10';

%market power: theta = 0.5
%this corresponds to Cournot competition with two high-quality plans
theta = .5;
out_baseline = equilibrium_function_v2(...
    v_sort,mc,q,theta,...
    bw,smooth_type,plot_switch,plot_number,plot_dir,plot_name,version); 

%% Equilibrium - Misperceptions

%sort matrix X in descending order by wtp and smooth
wtp_shift = 900;
X = [v_per + wtp_shift c_ins_mean_90-c_ins_mean_60];
X = sortrows(X,-1);

v_sort = smooth(X(:,1),bw,smooth_type);
mc = smooth(X(:,2),bw,smooth_type);
q = transpose((1:n_i)./n_i);

plot_switch = 1;
plot_number = 3;
plot_dir = '';
plot_name = 'all_rho_0';

theta = .5;
out_per = equilibrium_function_v2(...
    v_sort,mc,q,theta,...
    bw,smooth_type,plot_switch,plot_number,plot_dir,plot_name,version); 

%suplus under actual demand
X = [v_per v c_ins_mean_90-c_ins_mean_60];
X = sortrows(X, -1);

v_sort_perceived = X(:,1);
v_sort_actual = X(:,2);

p_star  = out_per(1,2);
q_star  = out_per(2,2);
costs   = out_per(3,2);
index = q_star*n_i;

producer_surplus            = p_star*index - costs;
consumer_surplus_perceived  = sum(v_sort_perceived(1:index)) - p_star*index;
consumer_surplus_actual     = sum(v_sort_actual(1:index)) - p_star*index;
total_surplus_perceived     = producer_surplus + consumer_surplus_perceived;
total_surplus_actual        = producer_surplus + consumer_surplus_actual;

fprintf('Producer surplus is %f \n',            producer_surplus);
fprintf('Perceived consumers surlus is %f \n',  consumer_surplus_perceived);
fprintf('Actual consumers surplus is %f \n',     consumer_surplus_actual);
fprintf('Perceived total surplus is %f \n',   total_surplus_perceived);
fprintf('Actual total surplus is %f \n',      total_surplus_actual);


